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Abstract — In this contribution, an algorithm for evaluating 
the capacity-achieving input covariance matrices for frequency 
selective Rayleigh MIMO channels is proposed. In contrast 
with the flat fading Rayleigh cases, no closed-form expressions 
for the eigenvectors of the optimum input covariance matrix 
are available. Classically, both the eigenvectors and eigenvalues 
are computed numerically and the corresponding optimization 
algorithms remain computationally very demanding. 

In this paper, it is proposed to optimize (w.r.t. the input covari- 
ance matrix) a large system approximation of the average mutual 
information derived by Moustakas and Simon. An algorithm 
based on an iterative water filling scheme is proposed, and its 
convergence is studied. Numerical simulation results show that, 
even for a moderate number of transmit and receive antennas, the 
new approach provides the same results as direct maximization 
approaches of the average mutual information. 

L Introduction 

When the channel state information is available at both 
the receiver and the transmitter of a MIMO system, the 
problem of designing the transmitter in order to maximize 
the (Gaussian) mutual information of the system has been 
addressed successfully in a number of papers. This problem is 
however more difficult when the transmitter has the knowledge 
of the statistical properties of the channel, a more realistic 
assumption in the context of mobile systems. In this case, 
the mutual information is replaced by the average mutual 
information (EMI), which, of course, is more complicated to 
optimize. 

The optimization problem of the EMI has been addressed 
extensively in the case of certain flat fading Rayleigh channels. 
In the context of the so-called Kronecker model, it has been 
shown by various authors (see e.g. |1| for a review) that 
the eigenvectors of the optimal input covariance matrix must 
coincide with the eigenvectors of the transmit correlation 
matrix. It is therefore sufficient to evaluate the eigenvalues of 
the optimal matrix, a problem which can be solved by using 
standard optimization algorithms. Similar results have been 
obtained for flat fading uncorrected Rician channels (121). 

In this paper, we consider this EMI maximization problem 
in the case of popular frequency selective MIMO channels 



(see e.g. f3l, f4l) with independent paths. In this context, the 
eigenvectors of the optimum transmit covariance matrix have 
no closed expressions, so that both the eigenvalues and the 
eigenvectors of the matrix have to be evaluated numerically. 
For this, it is possible to adapt the approach of [5j developed 
in the context of correlated Rician channels. However, the 
corresponding algorithms are computationally very demanding 
as they heavily rely on intensive Monte-Carlo simulations. 
We therefore propose to optimize the approximation of the 
EMI, derived by Moustakas and Simon (|4|), in principle valid 
when the number of transmit and receive antennas converge 
to infinity at the same rate, but accurate for realistic numbers 
of antennas. This will turn out to be a simpler problem. 
We mention that, while [4| contains some results related 
to the structure of the argument of the maximum of the 
EMI approximation, f4\ does not propose any optimization 
algorithm. 

We first review the results of |4| related to the large 
system approximation of the EMI. The expression of the 
approximation depends on the solutions of a non linear system. 
The existence and the uniqueness of the solutions is not 
addressed in |4|. As our optimization algorithm needs to solve 
this system, we clarify this crucial point. Next, we present our 
maximization algorithm of the EMI approximation. It is based 
on an iterative waterfiUing algorithm which, in some sense, 
can be seen as a generalization of |6l devoted to the Rayleigh 
context and of H devoted to the correlated Rician case: 
each iteration will be devoted to solve the above mentioned 
system of nonlinear equations as well as a standard waterfilling 
problem. It is proved that the algorithm converges towards the 
optimum input covariance matrix as long as it converge^. 

The paper is organized as follows. Section HI] is devoted 
to the presentation of the channel model, the underlying as- 
sumptions, the problem statement. The maximization problem 
of the EMI approximation is studied in section |lll] Numerical 
results are provided in section IIVI 

' Note however that we have been unable to prove formally its convergence. 



II. Problem statement 

A. General Notations 

In this paper, the notations s, x, M, stand for scalars, 
vectors and matrices, respectively. As usual, ||x|| represents 
the Euclidian norm of vector x, and ||M||, /o(M) and |M| 
respectively stand for the spectral norm, the spectral radius and 
the determinant of matrix M. The superscripts (.)-^ and (.)^ 
represent respectively the transpose and transpose conjugate. 
The trace of M is denoted by Tr(M). The mathematical 
expectation operator is denoted by E(-) 

All along this paper, r and t stand for the number of receive 
and transmit antennas. Certain quantities will be studied in the 
asymptotic regime t — >■ oo, r — >■ cx) in such a way thati/r — > 
c G (0,cx)). In order to simplify the notations, t ^ co should 
be understood from now on as i oo, ?- ^ cx) and t/r ^ 
c e (0, oo). 

Several variables used throughout this paper depend on 
various parameters, e.g. the number of antennas, the noise 
level, the covariance matrix of the transmitter, etc. In order to 
simplify the notations, we may not always mention all these 
dependencies. 

B. Channel model 

We consider a wireless MIMO link with t transmit and 
r receive antennas corrupted by a multi-paths propagation 
channel. The discrete-time propagation channel between the 
transmitter and the receiver is characterized by the input-output 
equation 

L 

y(n) = ^ H,s(n -/ + !) + n(n) = [H(z)]s(n) + n(n) (1) 

1=1 

where s(n) = (si(n), . . . , Sf (n))-^ represents the transmit 
vector at time n, y{n) = {yi{n), . . . ,yr{n))'^ the receive 
vector, and where n(n) is an additive Gaussian noise such that 
E(n(n)n(n)^) = cr^I. H(z) denotes the transfer function of 
the discrete-time equivalent channel defined by 
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(2) 



1=1 



Each coefficient H; is assumed to be a Gaussian random 
matrix given by 



(3) 



where W; is a r x t random matrix whose entries are inde- 
pendent and identically distributed complex circular Gaussian 
random variables, with zero mean and unit variance. The 
matrices C*^'' and C^'^ are positive definite, and account for 
the receive and transmit antenna correlation. We also assume 
that for each k ^ I, matrices and H; are independent. 

In the context of this paper, the channel matrices are as- 
sumed perfectly known at the receiver side. However, only the 
statistics of the (H.i)i=i....,l, i-S- matrices (C^'\ C'''')/=i_..._l, 
are available at the transmitter side. 



C. Ergodic capacity of the channel. 

Let Q(e^"'') be the t x i spectral density matrix of the 
transmit signal s(7i), which is assumed to verify the transmit 
power condition 

- f Tr(Q(e2--))dz. = 1 (4) 
^ Jo 

Then, the (Gaussian) ergodic mutual information /(Q(.)) 
between the transmitter and the receiver is defined as 

fi 
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I. + -.H(.)Q(.)H(.)^ 



dv 



(5) 



where Ew[-] = Ecw,)i=i £,[■]■ The ergodic capacity of the 
MIMO channel is equal to the maximum of J(Q(.)) over the 
set of all spectral density matrices satisfying the constraint 
dUl. The hypotheses formulated on the statistics of the channel 
allow however to limit the optimization to the set of positive 
matrices which are independent of the frequency v. This is 
because the probability distribution of matrix H(e^*'^'') is 
clearly independent of the frequency v. More precisely, the 
mutual information I(Q(.)) is also given by 



J(Q(.)) =Eh 



/'log I, + ^H(1)Q(.)H(1)^ 



dv 



where H = X^^Li H/ = H(l). Using the concavity of the 
logarithm, we obtain that 



/(Q(.)) < Eh 



log 



^H(l) 



Q{.)dv] H(l)^ 



We denote by C the cone of non negative hermitian matrices, 
and by Ci the subset of all matrices Q of C satisfying 
jTr(Q) = 1. If Q is an element of Ci, the mutual information 
/(Q) reduces to 



/(Q) = Eh 



log 



-HQH^ 



(6) 



It is strictly concave on the convex set Ci and reaches its 
maximum at a unique element Q, e Ci. It is clear that if 
Q^g2j7riy-j ^jjy spectral density satisfying (|4|i, then the matrix 
Q(e^*'^'^)dt^ is an element of Ci. Therefore, 

/(Q(.)) </(Q,) 

for each spectral density matrix verifying (|4|i. This shows that 
the maximum of function / over the set of all spectral densities 
satisfying (|4]i is reached on the set Ci. The ergodic capacity 
of the channel is thus equal to 



Ce — max /(Q) 

QeCi 



If the matrices fC^'^l 



coincide with a matrix C, matrix 
H follows a Kronecker model with transmit and receive 
covariance matrices X^f^i ^^'^ ^ respectively |8j. In this 
case, the eigenvectors of the optimum matrix coincide with 
the eigenvectors of J2i'=i The situation is similar if the 
transmit covariance matrices (C^'^)/=i^....l coincide. In the 
most general case, the eigenvectors of have however no 
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(7) 



closed form expression. The evaluation of Q* and of the chan- 
nel capacity Cb is thus a more difficult problem. A possible 
solution consists in adapating the Vu-Paulraj approach (ISJ) to 
the present context. However, the algorithm presented in 15] 
is very demanding since the evaluation of the gradient and the 
Hessian of /(Q) requires intensive Monte-Carlo simulations. 

D. The large system approximation of I{Q,). 

When t and r converge to oo while t/r c, c ^ (0, oo), ||4l 
showed that /(Q) can be approximated by /(Q) defined by (|7]l 
at the top of the page, where ((5i(Q), . . . , i5l(Q))"^ = <5(Q) 
and (^i(Q), . . . ,(5l(Q))^ — S{Q) are the positive solutions 
of the system of 2L equations: 

j Ki = !i{k) 

\ ki^ fi{K,Q) 

with K = (ki, . . . , kl)^ and k = {ki, . . . , ki)'^ , and 

/,(/i) = iTr[C(')T(/i)] 



(8) 



where 



/,KQ) = iTr QV2c(0Qi/2t(,c,Q) 



(9) 



(10) 



im is based on the replica method, a useful and simple trick 
whose mathematical relevance is not yet proved in the present 
context. However, using large random matrix technics similar 
to those of (191, H), it is possible to prove rigorously that, 
under mild technical extra assumptions, /(Q) — /(Q) + 0( j). 
This point is outside the scope of the present paper. 

We also mention that 14) assumed implicitely the exis- 
tence and the uniqueness of positive solutions of (O without 
justification. We therefore precise this important point, and 
also show that ((5;(Q))i=i,...,L and (5;(Q));=i,...,l can be 
evaluated using a fixed point algorithm. 

1) Existence: Using analytic continuation technique and 
results of [10], it can be shown that the following fixed point 
algorithm, initialized as follows, converges: 

. InitiaHzation: (jf"' > 0, 6^°^ > 0, I = 1, . . . , L. 



Evaluation of the ^["+^' and from (^("^ 

and S^"^ = Cs["\-.-,6''-Y-- 



(11) 



Besides, it can be proved that the limit of ((5^"-',^^ "*) when 
n — > 00 satisfies equation ([8]l, and that all the entries of this 
limit are positive. Hence, the convergence of the algorithm 
yields the existence of a solution to 



2) Uniqueness: In order to simplify the notations, we 
consider in this part the case Q = I. In order to address the 
general case, it is sufficient to change matrices (C('))i=i ... ^ 
into (Q^/2C(')Qi/2)i=i^...^i in what follows. Let (<5,^)'and 
{S' , 6 ) be two solutions of the canonical equation We 
denote (T, T) and (T',T') the associated matrices defined 
by ([Tol l. Introducing e — 6 — 6' = {ei, . . . , bl)'^ we have: 



ei 



iTr[c(')T(T' i_T-i)T' 



(12) 



k=l 



Similarly, with e = ^ — ^ = (ei, . . . , cl)^, 
^2 L 



efc 



t 



1=1 



J2iS'i-Si)Tr(c^'''>f&')f'] (13) 



And (O and (fTsT i can be written together as 



I cr2A(T,T') 
cr2A(f,f') I 



(14) 



with Aki{T,T) = iTr(C('=)TCWT') and Afe/(t,f') = 
iTr(C('=)tC(')f '). We will now prove that p{M) < 1, with 
M = cr4^(t,t')A(T,T'). This will imply that matrix M 
is invertible, and thus that e = e = 0. 



M.,| = 



4 



(15) 



Thanks to the inequality |Tr(AB)| < ^Tr(AAH)Tr(BBH) 
we have 



Tr(C('^')TC(^')f') < -y/Afcj(t,t)Afej(t',f') (16) 
Tr(C(^')TC(')T') 



< JA,,(T,T)A,,(T',T') (17) 



Using ( fT6] l and ( [Tt] ) in ( fTsT i gives 



|M,,| < v'Afe,(t)Afc,(f')A,,(T)A,,(T') 

i=i 

with A(f) A(f,T) and A(T) = A(T,T). And, using 
Cauchy-Schwarz inequality. 



\Mki\<a^ 



A,,(t)A,,(T) Afc, (f ')A,,(T') 



Hence, we have \M.ki\ < Pfei Vfc,Z, where the matrix P is 
defined by Pfe, = ^J (a4A(t)A(T))H ^ {a^A{f')A{T'))ki. 
Theorem 8.1.18 of flT\ then yields p{M) < p{P). Besides, 
Lemma 5.7.9 of [12] used on the definition of P gives: 



p(P) < (a4A(T)A(T)j (^a4A(T')A(T') 
We now introduce the following lemma: 

Lemma 1: p {a'^ k{T)A{T)^ < 1 
Proof: The 5i can be written as: 

6i = iTr(C(')TT-iT) 



t 



■Tr(C(')TT) 



— ^4Tr(C(')TC('=)T) 

fe=i 



" <5 









And ~5i = ^Tr(C(')tt) + ^ ^^^^ 4Tr(C(')tC('=)t) 
similarly, thus: 

A(T) 1 r 5 

A(t) \[5 

This equality is of the form u — Bu + v, where the entries 
of u and v are positive, and where the entries of B are non- 
negative. A direct application of Corollary 8.1.29 of ifTTI then 
implies p{B) < 1 - < 1. Noticing that (p(B))2 = 

p(cr4 A(f ) A(T)) ends the proof of Lemma [U □ 

We also have of course p{a^A{T')A{T')) < 1, so that 
finally: 

p{M) < p(P) < 1. 

III. Maximization algorithm 

Using the same methods as flSl section IV, the approxima- 
tion /(Q) can be shown to be a strictly concave function over 
the compact set Ci. Therefore it admits a unique argmax that 
we denote Q^. As Ci is convex, it is well known that is 
characterized by the property 



(V/(QJ,P-QJ <0 



(18) 



for each matrix P S Ci, where (V/(Q),P — Q) represents 
the limit of A-i (7(Q + A(P - Q)) - /(Q)) when A ^ 0, 
A > 0. We now consider the function V(Q, k, k) defined by: 



V(Q,k,k) =log|I + C(K)| +log|I + QC(K) 



(19) 



1=1 



where C(k) = J^Li and C{k) = Ya^i hC^^^ ■ Note 

that we have V(Q, <5(Q), ^(Q)) = T(Q). We have then the 
following result: 

Proposition 1: Denote d^, ~ ^(Q*) and S^, = ^(Q*). 
Matrix Q* is the solution of the standard waterfiUing problem: 
maximize over Q e Ci the function log |I + QC((5*)|. 

Proof: Due to lack of space, only the key points are given. We 
first remark that maximizing function Q i-^ log |I + QC(J»)| 



is equivalent to maximizing function Q i-^ V(Q,<5*,(5*) by 
(fT9] l. The proof is then based on the observation that 



dV 



(20) 
(21) 



are zero at point ((5(Q), <5(Q)). This implies that for each P e 
Ci, (V7(QJ,P-Q,> coincides with (VqV(Q^, (5*, ^,), P- 
Q^). As function Q V(Q, S^, 6^,) is strictly concave on Ci, 
(fTsT i implies that its argmax on Ci coincides with Q^. □ 

Proposition [T] shows that the optimum matrix is solution 
of a waterfilling problem associated to the covariance matrix 
C(^,). Although this result provides some insight on the 
structure of Q^, it cannot be used to evaluate it because matrix 
C(^*) depends itself of Q^. We now introduce an optimization 
algorithm of /(Q); the iterative scheme is the following: 

• Initialization: Qq = I 

• Evaluation of Qk fromQ^^i: {S^''\S ') is defined as 
the unique solution of ^ in which Q = Qk-i- Then 
Qfc is defined as the maximum of function Q i-^ 
log I + QCiS^''^) onCi. 

We now establish a result which shows that if the algorithm 
converges, then it converges towards Q,. 



Proposition 2: Assume that 



lim 6^''^ - S^''-^^ = lim S 



~Jk) ~(fc-l) 



(22) 



Then, the algorithm converges torwards matrix Q^. 

Proof: Due to the lack of space, we just outline the proof 
which is similar to the proof of Proposition 6 of [7|. As 
Ci is compact, we have just to verify that each convergent 
subsequence (Q^,(fc))fceN extracted from (Qfc)fcgN converges 
towards Q,. For this, we denote by Q^, ^ the limit of the above 
subsequence, and prove that this matrix verifies property ( fTSl ). 
We first remark that sequences ^'^('''+^ and S^'' ''^ converge 
towards vectors denoted S"^'* and ^'^ respectively. Moreover, 
{6'^'*,S ' ) is solution of system ([8]) in which matrix Q co- 
incides with ^. Therefore, using relations (|20] | and (l2Tl i as 
in the proof of Proposition [T] we obtain that (V/(Q^ P — 
Q^,*) coincides with (VqV(Q^^,, (5^,,,, J,^,,), P - Q^J. 
It remains to show that this term is negative for each P to 
complete the proof. For this, we use that Q.^^k) is the argmax 

over Ci of function Q ^ V(Q, d'''^''\d'^^''^). Therefore, 

(VQV(Q^i,(fc), ^^(fe), ^^(fe)), P - Q^(fe)) <0 (23) 

By (|22li, sequences {S^(^k))k>o and (^^(fc))fc>o converge to- 
wards S^'* and d respectively. Taking the limit of (1231 ) 
when fc oo shows that (Vq V(Q^ „, 5^,,,), P — 
Q^,») < as required. □ 




2 4 6 8 10 12 14 16 18 20 
SNR[dB] 



Fig. 1. Comparison with Vu-Paulraj algorithm 

To conclude, if the algorithm is convergent, that is, if the 
sequence of {Qk)keti converges towards a certain matrix, then 
the sl''^ — (5;(Qfc_i) and the (Jp-* — (5;(Qfc_i) converge as 
well when fc — > oo. ( |22] ) is verified, hence,if the algorithm 
is convergent it converges towards Q^. Although the con- 
vergence of the algorithm has not been proved, this result 
is encouraging and suggests that the algorithm is reliable. In 
particular, in all the conducted simulations the algorithm was 
converging. In any case, condition (l22l l can be easily checked. 
If it is not satisfied, it is possible to modify the initial point 
Qo as many times as needed to ensure the convergence. 

IV. Numerical Results 

We provide here some simulations results to visualize the 
impact of the transmit correlation optimization in a realistic 
context. We use the propagation model introduced in p|, in 
which each path corresponds to a scatterer cluster character- 
ized by a mean angle of departure and an angle spread. 

In the featured simulations, we consider a frequency se- 
lective MIMO system with i = r = 4, a carrier frequency 
of 2GHz, a number of paths L = 5. The paths share the 
same power, and their mean departure angles and angles 
spreads are given in Table U (in radians). In Figure [T] we 
have represented the true EMI /(I) (i.e. without optimization), 
and the optimized EMI /(Q,) (i.e. with an input covariance 
matrix maximizing the approximation). The EMI are evaluated 
by Monte-Carlo simulations, with 10^ channel realizations. 
The EMI optimized with Vu-Paulraj algorithm |5| is also 
represented for comparison. We fixed to 10 the number of 
iterations in [5:|. 

Figure [T] shows that maximizing /(Q) over the input 
covariance leads to significant improvement for /(Q). Our 
approach provides the same results as Vu-Paukaj's algorithm 
at high SNR, and performs even better elsewhere. Vu-Paulraj 's 
approach is penalized by the barrier method if the optimal 
input covariance is close to be singular, which is here the 
case if the SNR is not high enough. Moreover our algorithm 
is computationally much more efficient: in Vu-Paulraj 's algo- 
rithm the evaluation of the gradient and of the Hessian of /(Q) 
needs heavy Monte-Carlo simulations (10"' trials were used). 



TABLE I 

Paths angular parameters (in radians) 





I = 1 


I = 2 


I = 3 


I = 4 


I = 5 


mean departure angle 


6.15 


3.52 


4.04 


2.58 


2.66 


departure angle spread 


0.06 


0.09 


0.05 


0.05 


0.03 


mean arrival angle 


4.85 


3.48 


1.71 


5.31 


0.06 


arrival angle spread 


0.06 


0.08 


0.05 


0.02 


0.11 



TABLE II 
Average execution time (in seconds) 





L = 3 


L = 4 


L = 5 


Vu-Paulraj 


903 


1245 


1649 


New algorithm 


7,0.10-3 


7,4.10-3 


8,3.10-3 



Table gives for both algorithms the average execution time 
in seconds to obtain the input covariance matrix, on a 3.16GHz 
Intel Xeon CPU with 8GB of RAM, for a number of paths 
L = 3, L = 4 and L = 5. 

V. Conclusion 

In this paper we have addressed the evaluation of the 
capacity achieving covariance matrices of frequency selective 
MIMO channels. We have proposed to optimize a large system 
approximation of the EMI, and have introduced an attractive 
iterative algorithm. 
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